home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 08 - 1992 / 08.03 Jul 92 / Fast Random Numbers / PasRandomNumbers.p < prev    next >
Encoding:
Text File  |  1991-10-05  |  4.7 KB  |  164 lines  |  [TEXT/MPS ]

  1. UNIT RandomNumbers;
  2.  
  3. {  This unit implements a "minimal standard" random number
  4.    generator using 32-bit integer arithmetic.  This is not
  5.    the "best" random number generator possible, but it is
  6.    simple, quick, and produces numbers which are "random"
  7.    enough for most purposes.
  8.  
  9.    It is based on the "parametric multiplicative linear con-
  10.    gruential algorithm" which was originally proposed by D.
  11.    H. Lehmer in 1951.  This algorithm generates a sequence
  12.    of integers z[1], z[2], z[3]... by repeatedly carrying 
  13.    out the following calculation:
  14.  
  15.                 z[n+1] = (a * z[n]) mod m
  16.  
  17.    where "a" and "m" are suitable constants.  The successive
  18.    "z"s are stored in the variable "seed" in this unit.  
  19.    They all lie in the range [0..m-1], so dividing them by
  20.    "m" gives real numbers which lie in the range [0.0..1.0].
  21.  
  22.    In 1969, Lewis, Goodman and Miller proposed the choice of
  23.    a = 16807 and m = 2**31 - 1 = 2147483647 as particularly
  24.    suited for a machine with a 32-bit word length (in their
  25.    case, the IBM System/360).  The main problem is that the
  26.    product a * z[n] can be more than 32 bits long, leading
  27.    to integer overflow on many systems (including the Mac).
  28.  
  29.    In 1979, G. L. Schrage devised an implementation of the
  30.    LG&M method which is guaranteed not to produce integer
  31.    overflow on any system with a 32-bit word length.  It is
  32.    this implementation which is used here, in procedure
  33.    "UpdateSeed".
  34.  
  35.    Source:  Stephen K. Park and Keith W. Miller, "Random
  36.    Number Generators:  Good Ones Are Hard to Find", Communi-
  37.    cations of the ACM, vol. 31, p. 1192  (October 1988).
  38.    
  39.    This unit was written in MPW Pascal, v3.2.
  40.  
  41.                            Jon Bell
  42.               Dept. of Physics & Computer Science
  43.                      Presbyterian College
  44.                        Clinton SC 29325
  45.                     CompuServe: #70441,353
  46.                           October 1991                     }
  47.  
  48. INTERFACE
  49.  
  50. {----------------------------------------------------------}
  51.  
  52. procedure InitRandomSeed (newSeed : longint);
  53.  
  54. {  Initializes the random number seed to "newSeed".  You
  55.    must call this procedure once, at the beginning of your
  56.    program, before you use any of the following functions.
  57.    As far as randomness is concerned, it makes no difference
  58.    what value you use for "newSeed", so long as it isn't
  59.    zero.  Using different seeds merely gives you different
  60.    sequences of "random" numbers.  Using the same seed each
  61.    time you run the program gives you the same sequence of
  62.    "random" numbers each time, which may be useful for
  63.    debugging. }
  64.  
  65. {----------------------------------------------------------}
  66.  
  67. function RandomSeed : longint;
  68.  
  69. {   Returns the current value of the random number seed.   }
  70.  
  71. {----------------------------------------------------------}
  72.  
  73. function RandomReal : extended;
  74.  
  75. {  Returns a real number, uniformly "randomly" selected 
  76.    from the interval (0.0,1.0).  Note that this is an
  77.    "open" interval, that is, it does _not_ include 0.0 or
  78.    1.0.                                                    }
  79.  
  80. {----------------------------------------------------------}
  81.  
  82. function RandomInteger (max : longint) : longint;
  83.  
  84. {  Returns an integer, uniformly "randomly" selected from
  85.    the interval [1,max].  Note that this is a "closed"
  86.    interval, that is, it _does_ include 1 and max.         }
  87.  
  88. {----------------------------------------------------------}
  89. {----------------------------------------------------------}
  90.  
  91. IMPLEMENTATION
  92.  
  93. const
  94. a = 16807;
  95. m = 2147483647;
  96. q = 127773;       { m div a }
  97. r = 2836;         { m mod a }
  98.  
  99. {  NOTE:  Other possible values for these constants, which
  100.    may give "better" results, are:
  101.  
  102.       a = 48271,  q = 44488,  r = 3399  or
  103.       a = 69621,  q = 30845,  r = 23902,
  104.  
  105.    keeping m = 2147483647.                                 }
  106.  
  107. var
  108. seed : longint;
  109.  
  110. {----------------------------------------------------------}
  111.  
  112. procedure InitRandomSeed (newSeed : longint);
  113.  
  114. begin
  115. seed := newSeed
  116. end;
  117.  
  118. {----------------------------------------------------------}
  119.  
  120. function RandomSeed : longint;
  121.  
  122. begin
  123. RandomSeed := seed
  124. end;
  125.  
  126. {----------------------------------------------------------}
  127.  
  128. {private} procedure UpdateSeed;
  129.  
  130. var
  131. lo, hi, test : longint;
  132.  
  133. begin
  134. hi := seed div q;
  135. lo := seed mod q;
  136. test := a * lo - r * hi;
  137. if test > 0 then
  138.     seed := test
  139. else
  140.     seed := test + m
  141. end;
  142.  
  143. {----------------------------------------------------------}
  144.  
  145. function RandomReal : extended;
  146.  
  147. begin
  148. UpdateSeed;
  149. RandomReal := seed / m
  150. end;
  151.  
  152. {----------------------------------------------------------}
  153.  
  154. function RandomInteger (max : longint) : longint;
  155.  
  156. begin
  157. UpdateSeed;
  158. RandomInteger := (seed mod max) + 1
  159. end;
  160.  
  161. {----------------------------------------------------------}
  162.  
  163. END.
  164.